from PIL import Image, ImageOps, ImageDraw
import numpy as np
from matplotlib import pyplot as plt
import colorsys
import math
import pandas as pd

Data undestanding and pre - processing

input_filename_1 = "population-density-map.bmp"
input_filename_2 = "elevation1x1_new-mer-bleue.bmp"
im_popul = Image.open(input_filename_1)
im_elev = Image.open(input_filename_2)

display(im_popul)

display(im_elev)

im_elev.size , im_popul.size #size of each imagen

((4901, 4251), (4830, 3510))

Population map

width=im_popul.size[0]
heigth=im_popul.size[1]
colors = im_popul.getcolors(width*heigth)
print('Nb of different colors: %d' % len(colors))

Nb of different colors: 32

According to the information provided, the population density file has a red pixel for the center of Rize and a green pixel for the center of Brest, all other pixels are black and white. Maximum density (absolute white color) is 3000 inhabitants per km2 (every pixel is a square kilometer).

def hexencode(rgb):
    r=rgb[0]
    g=rgb[1]
    b=rgb[2]
    return '#%02x%02x%02x' % (r,g,b)

for idx, c in enumerate(colors):
    plt.bar(idx, c[0], color=hexencode(c[1]),edgecolor=hexencode(c[1]))


plt.show()

colors

[(91189, (255, 255, 255)),
 (28235, (246, 246, 246)),
 (36939, (238, 238, 238)),
 (47019, (230, 230, 230)),
 (67539, (222, 222, 222)),
 (74903, (213, 213, 213)),
 (87366, (205, 205, 205)),
 (99977, (197, 197, 197)),
 (113816, (189, 189, 189)),
 (147283, (180, 180, 180)),
 (151062, (172, 172, 172)),
 (169603, (164, 164, 164)),
 (183049, (156, 156, 156)),
 (218401, (148, 148, 148)),
 (194375, (139, 139, 139)),
 (195366, (131, 131, 131)),
 (177610, (123, 123, 123)),
 (157859, (115, 115, 115)),
 (151124, (106, 106, 106)),
 (101929, (98, 98, 98)),
 (87459, (90, 90, 90)),
 (49058, (82, 82, 82)),
 (72353, (74, 74, 74)),
 (21953, (65, 65, 65)),
 (7757, (57, 57, 57)),
 (3000, (49, 49, 49)),
 (1357, (41, 41, 41)),
 (301, (32, 32, 32)),
 (12252255, (24, 24, 24)),
 (1963161, (16, 16, 16)),
 (1, (0, 255, 0)),
 (1, (255, 0, 0))]

There are 91189 pixel of complete white color, in other words where there exist 91189* 3000 inhabitants that can help to prevent fight

pixels_pob = np.array(im_popul) 

print(pixels_pob.shape) #3510 is the height and 4830 the width

(3510, 4830, 3)

# With the image im, let's generate a numpy array to manipulate pixels
p = np.array(im_popul) 
print(p.shape)

(3510, 4830, 3)

mask = np.all(p == (0, 255, 0), axis=-1)
z = np.transpose(np.where(mask))
print("Coordinates (x,y) of the green pixel (Rinze): (%d,%d)" % (z[0][1],z[0][0]))
mask = np.all(p == (255, 0, 0), axis=-1)
z = np.transpose(np.where(mask))
print("Coordinates (x,y) of the red pixel (Brest): (%d,%d)" % (z[0][1],z[0][0]))

Coordinates (x,y) of the green pixel (Rinze): (4426,2108)
Coordinates (x,y) of the red pixel (Brest): (669,1306)

grayim_pop = im_popul.convert("L")
display(grayim_pop)

colors = grayim_pop.getcolors(width*heigth)
print('Nb of different colors: %d' % len(colors))
# With the image im, let's generate a numpy array to manipulate pixels
pop = np.array(grayim_pop) 
# plot the histogram. We still have a lot of dark colors. Just to check ;-)
plt.hist(pop.ravel())

Nb of different colors: 32

(array([1.4215717e+07, 1.2114000e+04, 1.4336500e+05, 3.4051200e+05,
        5.3083500e+05, 5.9582600e+05, 4.6794800e+05, 3.0115900e+05,
        1.8946100e+05, 1.5636300e+05]),
 array([ 16. ,  39.9,  63.8,  87.7, 111.6, 135.5, 159.4, 183.3, 207.2,
        231.1, 255. ]),
 <BarContainer object of 10 artists>)

# from gray colors to density
density = pop/255.0
# plot the histogram. We still have a lot of dark colors. Just to check ;-)
plt.hist(density.ravel())

(array([1.4215717e+07, 1.2114000e+04, 1.4336500e+05, 3.4051200e+05,
        5.3083500e+05, 5.9582600e+05, 4.6794800e+05, 3.0115900e+05,
        1.8946100e+05, 1.5636300e+05]),
 array([0.0627451 , 0.15647059, 0.25019608, 0.34392157, 0.43764706,
        0.53137255, 0.62509804, 0.71882353, 0.81254902, 0.90627451,
        1.        ]),
 <BarContainer object of 10 artists>)

print("Population density in brest: {}".format(density[1306,669]))
print("Population density in Rinze: {}".format(density[2108,4426]))#Poblacional density in brest

Population density in brest: 0.2980392156862745
Population density in Rinze: 0.5882352941176471

Elevation map

width=im_elev.size[0]
heigth=im_elev.size[1]
colors = im_elev.getcolors(width*heigth)
print('Nb of different colors: %d' % len(colors))
# To plot an histogram
from matplotlib import pyplot as plt
def hexencode(rgb):
    r=rgb[0]
    g=rgb[1]
    b=rgb[2]
    return '#%02x%02x%02x' % (r,g,b)

for idx, c in enumerate(colors):
    plt.bar(idx, c[0], color=hexencode(c[1]))

plt.show()

Nb of different colors: 189

The blue corresponds to the sea, being this the distance 0m and the darker brown corresponds to the highest point, which corresponds to the mont blanc located at approximately 4,809 m, it is the highest point in Europe.

#HSV - for hue, saturation, value

def HSVColor(img):
    
    r,g,b = img.split()
    Hdat = []
    Sdat = []
    Vdat = [] 
    for rd,gn,bl in zip(r.getdata(),g.getdata(),b.getdata()) :
        h,s,v = colorsys.rgb_to_hsv(rd/255.,gn/255.,bl/255.)
        Hdat.append(int(h*255.))
        Sdat.append(int(s*255.))
        Vdat.append(int(v*255.))
    r.putdata(Hdat)
    g.putdata(Sdat)
    b.putdata(Vdat)
    return Image.merge('RGB',(r,g,b))

Taking the Hue saturation, 0 red and 240 blue we can obtain another representation of the image

image_hsv = HSVColor(im_elev)

pixels_Hue = np.array(image_hsv) 
plt.imshow(pixels_Hue[:,:,0], cmap ="gray")

<matplotlib.image.AxesImage at 0x7f66d46b2fa0>

pix_H = pixels_Hue[:,:,0] #Dimension Hue
unique= np.unique(pix_H)
print(len(unique))
unique

75

array([  0,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,  14,  15,
        16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,  28,
        29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
        42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,
        55,  56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,
        68,  69,  70,  71,  72,  73,  74,  75,  76, 138], dtype=uint8)

x_1,y_1 = 0 , 255 # 0 :lower height and 255:highest altitud
x_2,y_2 = 138 , 0 # 138 : Blue see color 0: red mountain color

m = (y_2 - y_1) / (x_2 -x_1)
print(m)
aux_x1 = np.full((pix_H.shape[0], pix_H.shape[1]), x_1)
aux_y1 = np.full((pix_H.shape[0], pix_H.shape[1]), y_1)
pix_final = m * (pix_H - aux_x1 ) + aux_y1

-1.8478260869565217

plt.imshow(pix_final, cmap ="gray")

<matplotlib.image.AxesImage at 0x7f66d0696580>

im = Image.fromarray(pix_final)
im = im.convert('L') 
im.save("Elevation_bw_2.bmp") #Export to use matlabcode

input_filename_3 = "elevation_fixed.bmp"
im_elevT1 = Image.open(input_filename_3)
#im_elevT1 = im_elev
width=im_elevT1.size[0]
heigth=im_elevT1.size[1]
display(im_elevT1)

colors= im_elevT1.getcolors(width*heigth)
print('Nb of different colors: %d' % len(colors))

Nb of different colors: 256

im_elevT1.size , im_popul.size #size of each imagen

((4830, 3510), (4830, 3510))

pop = np.array(grayim_pop) 
pop.shape

(3510, 4830)

elevation = np.array(im_elevT1)
elevation.shape

(3510, 4830)

plt.imshow(pop[:3000,:], cmap ="gray")

<matplotlib.image.AxesImage at 0x7f66d48c8640>

plt.imshow(elevation[:3000,:], cmap ="gray")

<matplotlib.image.AxesImage at 0x7f66d471d0a0>

Pop_final = pop[:3000,:]
Ele_final = elevation[:3000,:]

pop_density = Pop_final /255 # 1 = 3000 inhabitans 
ele_density = Ele_final /255 # 1 = 4809 m altitud

plt.hist(pop_density.ravel())

(array([1.1891935e+07, 1.1257000e+04, 9.7343000e+04, 3.2558700e+05,
        5.2126100e+05, 5.8097700e+05, 4.5177800e+05, 2.8796500e+05,
        1.7934500e+05, 1.4255200e+05]),
 array([0.0627451 , 0.15647059, 0.25019608, 0.34392157, 0.43764706,
        0.53137255, 0.62509804, 0.71882353, 0.81254902, 0.90627451,
        1.        ]),
 <BarContainer object of 10 artists>)

plt.hist(ele_density.ravel())

(array([5300111.,   19073.,   18201.,   18559.,  925879., 2305608.,
        1751305., 1434374., 1766935.,  949955.]),
 array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
 <BarContainer object of 10 artists>)

print("Population density in brest: {}".format(pop_density[1306,669]))
print("Population density in Rinze: {}".format(pop_density[2108,4426]))

Population density in brest: 0.2980392156862745
Population density in Rinze: 0.5882352941176471

print("Altitud density in brest: {}".format(ele_density[1306,669]))
print("Altitud density in Rinze: {}".format(ele_density[2113,4426]))

Altitud density in brest: 0.5058823529411764
Altitud density in Rinze: 0.592156862745098

a = np.array([[1,1],[1,1]])
b = np.array([[2,1/2],[1/2,1/2]])

a*b

array([[2. , 0.5],
       [0.5, 0.5]])

Graph construction

pop_density.max()

1.0

from scipy import signal

kernel_1 = np.ones((15,15)) / (15*15)
conv_pop = signal.convolve2d(pop_density, kernel_1 , mode='valid')
conv_ele = signal.convolve2d(ele_density, kernel_1 , mode='valid')

plt.imshow(np.absolute(grad), cmap='gray')

<matplotlib.image.AxesImage at 0x7f6636a4c5b0>

conv_pop.shape , pop_density.shape , conv_ele.shape , ele_density.shape

((2986, 4816), (3000, 4830), (2986, 4816), (3000, 4830))

((pop_density.shape[0] - kernel_1.shape[0] )// 15)+1  , ((pop_density.shape[1] - kernel_1.shape[1] )// 15)+1

(200, 322)

conv_pop_final=conv_pop[list(range(0,conv_pop.shape[0],15)),:]
conv_pop_final=conv_pop_final[:,list(range(0,conv_pop.shape[1],15))]

conv_ele_final=conv_ele[list(range(0,conv_ele.shape[0],15)),:]
conv_ele_final=conv_ele_final[:,list(range(0,conv_ele.shape[1],15))]

conv_pop_final.shape , conv_ele_final.shape 

((200, 322), (200, 322))

conv_pop_final.min() 

0.08658823529411776

np.count_nonzero(conv_pop_final < 0.09)

14992

np.count_nonzero(conv_ele_final == 0.0)

21068

conv_pop_final[:200,100:200]

array([[0.08937691, 0.08923747, 0.32095861, ..., 0.0962963 , 0.10086275,
        0.08909804],
       [0.11105882, 0.1292549 , 0.37233987, ..., 0.09021351, 0.08854031,
        0.08965577],
       [0.09221786, 0.08951634, 0.21605229, ..., 0.08965577, 0.09118954,
        0.09160784],
       ...,
       [0.08979521, 0.0964183 , 0.08979521, ..., 0.09049237, 0.08965577,
        0.08923747],
       [0.08854031, 0.08909804, 0.08965577, ..., 0.09035294, 0.09063181,
        0.08923747],
       [0.09091068, 0.08909804, 0.09035294, ..., 0.08951634, 0.08951634,
        0.08993464]])

[1306 //15,669//15] ,[2108//15,4426//15] #Brest and Rinze

([87, 44], [140, 295])

conv_pop_final[87,44] , conv_pop_final[140,295]

(0.5107276688453152, 0.3797298474945531)

conv_pop_final.min() , conv_pop_final.max()

(0.08658823529411776, 0.9982222222222237)

conv_ele_final.min(), conv_ele_final.max()

(0.0, 0.9978736383442283)

def cal_ang (x,y,alt_x,alt_y):
    
    alpha = 0
    d = math.sqrt(pow(x[0] - y[0],2) + pow(x[1] - y[1],2))
    #print(d)
    h = alt_y - alt_x
    #print(h)
    ang = math.degrees(math.atan(h/d))
    #print(ang)
    if ang <= 0:
        alpha = 1
    elif (ang > 0) and (ang < 10):
        alpha = 1 - (ang/10)
    else :
        alpha = 0 
        
    return alpha

x = (0,0)
y = (0,1)
alt_x = 0.9
alt_y = 1
ang = cal_ang (x,y,alt_x,alt_y)
print(ang)

1.0
0.09999999999999998
5.710593137499642
0.4289406862500358

Data = pd.DataFrame(columns = ['node_id','x','y','edge_attr','human_pop','zombie_pop'])

Data= Data.append({'node_id': 1,'x':0,'y':0,'edge_attr':0.1,'human_pop':100,'zombie_pop':1}, ignore_index=True)

Data= Data.append({'node_id': 2,'x':0,'y':0,'edge_attr':0.1,'human_pop':100,'zombie_pop':1}, ignore_index=True)

Data

	node_id	x	y	edge_attr	human_pop	zombie_pop
0	1.0	0.0	0.0	0.1	100.0	1.0
1	2.0	0.0	0.0	0.1	100.0	1.0

import networkx as nx

nodes_x = 200
nodes_y = 322

G = nx.grid_2d_graph(nodes_x,nodes_y).to_directed()
new_edges = [[((i,j),(i-1,j-1)),((i,j),(i-1,j+1)),((i,j),(i+1,j-1)),((i,j),(i+1,j+1))] 
                for i in range(0,nodes_x) for j in range(0,nodes_y)]
#new_edges += [[((1,0),(0,1)), ((0,nodes-2),(1,nodes-1)), ((nodes-2,0),(nodes-1,1)), ((nodes-2,nodes-1),(nodes-1,nodes-2))]]
G.add_edges_from(sum(new_edges,[]))

G_aux = G.copy()
#### To fix diagonals

for node in G.nodes():
    if (node[0] < 0) or (node[1] < 0 ) or (node[0] == nodes_x) or (node[1] == nodes_y ) :
        G_aux.remove_node(node)

plt.figure(figsize=(6,6))
pos = {(x,y):(y,-x) for x,y in G_aux.nodes()}
# nx.draw(G_aux, pos=pos, 
#         node_color='lightgreen', 
#         with_labels=False,
#         node_size=50)

nx.draw(G_aux, pos=pos, node_size=10)

# for node in G_aux.nodes():
#     print(node)
#     print(conv_pop_final[node[0],node[1]]*3000)
#     print(conv_ele_final[node[0],node[1]])


nx.set_node_attributes(G_aux, {node: {'node_id': 'U' + str(i), 'human_pop': int(conv_pop_final[node[0],node[1]]*3000), 'zombie_pop': 0} 
                                for i, node in enumerate(G_aux.nodes)})

#Node_id=nx.get_node_attributes(G_aux,'node_id')
# for i, n in enumerate(G_aux.nodes):
#     print(i)
#     print(n)

[n for n in G_aux.neighbors((1,1))]

[(0, 1), (2, 1), (1, 0), (1, 2), (0, 0), (0, 2), (2, 0), (2, 2)]

dic_edges= {e:cal_ang(e[0],e[1],conv_ele_final[e[0]],conv_ele_final[e[1]]) for e in G_aux.edges}

nx.set_edge_attributes(G_aux, dic_edges, name = 'elev_factor')

elev=nx.get_edge_attributes(G_aux,'elev_factor')

See_location= np.transpose(np.nonzero(conv_ele_final ==0))
See_l = [(g[0],g[1]) for g in See_location]

G_aux_2 = G_aux.copy()

for node in G_aux.nodes():
    if node in See_l:
        G_aux_2.remove_node(node)

plt.figure(figsize=(6,6))
pos = {(x,y):(y,-x) for x,y in G_aux_2.nodes()}
# nx.draw(G_aux_2, pos=pos, 
#         node_color='lightgreen', 
#         with_labels=True,
#         node_size=600)
nx.draw(G_aux_2, pos=pos, 
        node_color='lightgreen',
        node_size=10 , width=0.5)

loc_brest= (87,44)
loc_Rinze = (140, 295)

nx.get_node_attributes(G_aux_2, 'human_pop')[loc_Rinze]

1139

str(loc_Rinze)

'(140, 295)'

nx.set_node_attributes(G_aux_2, {loc_Rinze: int(conv_pop_final[loc_Rinze[0],loc_Rinze[1]]*3000)}, name = 'zombie_pop')
nx.set_node_attributes(G_aux_2, {loc_Rinze: 0}, name = 'human_pop')

nx.get_node_attributes(G_aux_2, 'node_id')[loc_Rinze], nx.get_node_attributes(G_aux_2, 'zombie_pop')[loc_Rinze] , nx.get_node_attributes(G_aux_2, 'human_pop')[loc_Rinze]

('U45375', 1139, 0)

nx.get_node_attributes(G_aux_2, 'node_id')[loc_brest], nx.get_node_attributes(G_aux_2, 'zombie_pop')[loc_brest] , nx.get_node_attributes(G_aux_2, 'human_pop')[loc_brest]

('U28058', 0, 1532)

#nx.write_gexf(G_aux, "All_nodes_graph.gexf")
nx.write_gexf(G_aux_2, "without_see_nodes_graph_ini_2.gexf")

import pickle

a_file = open("Positions_without_see_graph_2.pkl", "wb")
pickle.dump(pos, a_file)
a_file.close()

len(G_aux_2.nodes()) , len(G_aux_2.edges())

(43332, 335912)

G_zombie = nx.read_gexf("without_see_nodes_graph_ini_2.gexf")

len(G_zombie.nodes()) , len(G_zombie.edges())

(43332, 335912)

elev_2=nx.get_edge_attributes(G_zombie,'elev_factor')